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Spontaneous symmetry breaking is well understood through the classical “Mexican Hat” picture, 
which describe many quantum phases of matter. Recently, several new classes of quantum phases of 
matter, such as topological orders and symmetry protected topological (SPT) orders, were discov¬ 
ered. In an attempt to address the transitions between all those phases of quantum matter under 
the same framework, we introduced an analogous yet very simple picture for phase transitions in 
the context of tensor-networks. Using a very simple iteration process, we found that both symmetry 
breaking and some topological phase transitions (for topological orders described by gauge theory 
and ID SPT orders) could be marked by a sudden change in the symmetry structure of the so-called 
“environment matrix”. In this process, the environment matrix serves as an “order parameter” that 
captures patterns of entanglement in topological phases. The symmetry change in the environment 
matrix is very much like the symmetry breaking of conventional order parameters. We applied this 
method to both the transverse Ising model (ID and 2D honeycomb), spin-1 model (ID), and the 
Toric Code model in a magnetic field (2D honeycomb), and explored the corresponding symmetry 
structure changes in their environment matrices in details. With just a few variational parameters 
and a few minutes’ run time on a laptop, we could get the corresponding phase transition points 
within a few percent error compared with the Quantum Monte Carlo results. 


I. INTRODUCTION 

In recent years, with the discoveries of quantum Hall 
states^’^ and topological insulators,the field of con¬ 
densed matter physics is focusing more and more on 
topological phases of matter. Lots of progress has been 
made in the classification of topological order'^'^^ in in¬ 
teracting bosonic/fermionic systems through tensor net¬ 
work representation of many-body wave function and 
the associated fixed-point tensors under wave function 
renormalization,which lead to tensor category the¬ 
ory of topological order. 18 presence of sym¬ 

metry, tensor network and group cohomologyalso 
lead to a classification of symmetry protected topologi¬ 
cal (SPT) order. 

However, an important question is how to deter¬ 
mine the topological order or SPT order carried by a 
generic wave function or a generic tensor network wave 
function. ^bi2,14,23-32 tensors in different generic ten¬ 

sor network can look similar, but represent different topo- 
logical/SPT prders. This is because topological order 
is highly non-local. All its features, including ground- 
state degeneracy, braidings and statistics of the quasi¬ 
particles, topological entanglement entropy are global 
features. One can only see those features after perform¬ 
ing the wave function renormalization. This makes tra¬ 
ditional theory of using local “order parameters” to de¬ 
scribe topological/SPT orders impossible. 

Another difficulty to read topological/SPT orders from 
the local tensor is that although the existing matrix- 
product representation has reached great success in 
ID, its higher-dimension extension is still a numerically 
formidable task, and many 2D tensor-network renor¬ 


malization scheme face the infamous “corner double¬ 
line” problem,which tensor-network renormalization 
quickly break down after a few iterations. Thus a compu¬ 
tationally efficient tensor-network method to implement 
RG is badly desired. It was in view of this that we de¬ 
veloped our “mean-field” approach based on the environ¬ 
ment matrix. 

This paper is structured as follows: In section H, we 
first introduce the concept of “environment matrix” and 
outline how this method is applied in ID, with the ex¬ 
ample of the transverse Ising model. In section HI, we 
make some detailed emphasis on the symmetry structure 
of the environment matrix, which leads to a characteri¬ 
zation of different phases. In section IV and V, we detail 
how to detect different phases without knowing the sym¬ 
metry structure, which makes our method immediately 
applicable to existing ID numerical methods in identify¬ 
ing different SPT phases. Finally in section VI and VH, 
we generalize this method to 2D, and apply it to both the 
transverse Ising model and the Toric-Gode model with a 
B-field on a honeycomb lattice. 


II. ID ENVIRONMENT TENSOR METHOD 

The environment tensor method has been widely ap¬ 
plied in ID systems through the study of Matrix-Product 
States (MPS).’*^’^® Consider the ID transverse Ising 
model on an infinite lattice: 
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FIG. 1. A matrix-product state. All the physical sites are 
represented by dots, and physical/internal degrees of freedom 
by vertical/horizontal lines. 





FIG. 3. A self-consistent condition that environment matrix 
E must satisfy, where A is a scaling factor. The part enclosed 
by the shaded area is the so-called double-tensor Tbp,aa- 






Ev 


FIG. 2. Approximate the energy through the environment 
matrix, E. Note that because the system is translation invari¬ 
ant, we only need to calculate H for two neighboring sites. In 
ID, by assuming M to be left-right symmetric, the environ¬ 
ment matrices on the left and right will be the same. 


where ct’s are the regular Pauli matrices. Recall that the 
wave function of a ID system could always be written 
into a matrix-product form; in particular, if the system 
is translationally symmetric, then we have (in Figure 1): 

= E E (2) 

{mi} {at} J 

where matrices M’s are independent of the site labels i 
and are labeled by the physical degrees of freedoms m^. 

In ID, the environment tensor method is essentially 
a variational calculation based on the above matrix- 
product state, with matrices M’s as variational param¬ 
eters. We use the matrices M’s to obtain the average 
energy (see the top of Fig. 2). We minimize the the 
average energy to obtain M’s. 

The calculation of average energy is actually a finite 
calculation. The key is to use “environment matrix” to 
capture the contributions from far-away sites. This is 
graphically shown in the bottom of Fig. 2. As can be 
seen, there are two environment matrices, one on each 
side of the ID chain. 

So in the actual environment tensor method, we use 
the matrices M’s to obtain the environment matrices E, 
and then use the matrices M’s and environment matrices 
E’s to obtain the average energy. We then minimize the 
average energy to obtain M’s (and E's). 

For fixed matrices M’s, the environment matrix could 
be obtained through iterations. As shown in Figure 3, 
starting from some random initial values Eq that satis¬ 
fies TtEqEq = 1, we can update the environment matrix 
using a “double-tensor”, which is formed by two M ma¬ 
trices with physical indices contracted. After applying 
the “double-tensor”, Eq is changed to XiEi where Ei 
satisfies Tri?|i?i = 1 and Ai is a scaling factor. After 
iterating enough number of times, a final stable “envi¬ 


ronment matrix” Eao = E and a final stable scaling fac¬ 
tor Aoo = A would be reached. Note that this process, 
after viewing the environment as a vector and the double¬ 
tensor as an operator, is essentially equivalent to picking 
out the eigenvector with the largest absolute value of the 
eigenvalues of the double tensor. In this way, for each M, 
we can obtain the corresponding environment matrix E 
through iterations, and by applying E both on the left 
and on the right (see Figure 2), we can get the total en¬ 
ergy. The variational calculation could then be carried 
out for different values oih/ J, and a phase diagram could 
then be obtained. 

More specifically, we require our matrix-product state 
to have a Z 2 symmetry that corresponds to spin up-down 
flipping, for both the symmetry-breaking and the sym¬ 
metric phases. So even in symmetry breaking phase, we 
choose the ground state to be, say, (|tt ...) + ...)) 

when h/ J —>■ 0. Note that this is different from tradi¬ 
tional symmetry breaking description, where the ground 
state spontaneously picks one of the ferromagnetic states. 

Recall that on-site symmetry of the ground state re¬ 
quires matrices M’s to transform in a special way w.r.t. 
symmetry;^® 


^ (3) 

m' 

Here, m is the spin index, matrix g^m' represents the 
on-site symmetry and acts on the spin basis, 0 is a phase 
factor (set to 0 in this paper), Ug is a unitary matrix act¬ 
ing on internal degrees of freedom, and forms a projective 
representation of the symmetry group 

For internal dimension D (dimension of M) being 2, 

we can choose Ug = , then equation (3) reduces 

to 


M^ = UlM^Ug, 


( 4 ) 


and we thus have: M^ = 


I ^ and M4 = 


-b 


Here M’s are symmetric because of left-right symmetry, 
and a, b, c are free variational parameters. 

With the above symmetry analysis in mind, numerical 
simulation could be run on our ID Ising model. Following 
the previous discussion, for each h/ J G [0, 00 ] in equation 
(1), we minimize the energy by varying M’s satisfying 
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FIG. 4. Energy as a function of h/J. The h- and J- terms are 
also individually plotted in the graph, so the phase transition 
could be easily spotted at h/J = 0.83. This result is obtained 
when internal dimension D=2, and we’ve chosen the grid so 
that sample points are denser close to the transition point. 
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FIG. 5. Energy plot when internal dimension increased to 
D=4. We can see the phase transition point occured around 
h/J = 0.97, a big improvement from Fig. 4. 


equation (4). By plotting the two energy terms, a phase 
diagram is obtained (see Fig. 4). For internal dimension 
D = 2, the phase transition occurred at h/J = 0.83, with 
an error of 17%. 


We could easily improve the result by increasing the 
internal dimension. For internal dimension Z? = 4, we 

can choose Ug in (4) to be | ^ j, where I is the 2x2 


Identity matrix. The most general symmetric M satisfy¬ 
ing (4) has 10 variational parameters. Following the same 
variational procedure, we can get the energy plot shown 
in Fig. 5. The phase transition occurred at h/J = 0.97, 
with a mere 3% error. 


Note that in both calculations, we used symmetric ma¬ 
trices M’s with all real parameters. The typical runtime 


on a laptop was just a few seconds in both cases. 


III. SYMMETRY STRUCTURE OF THE 
ENVIRONMENT MATRIX 


From the above plot, we see that there is a phase tran¬ 
sition at h/J « 0.83. To understand the phases on the 
two sides of the transition, let us choose another basis 


2 \^ a — c a + c — 2bJ 
= WM^W^ = - ( 5 ) 

2 \ a - c a + c + 2bj ^ ^ 


where W = 2 



In the new basis the meaning 


of the M’s is more clear. 

When 6 = 0, M^ = M'^, and the MPS is a pure product 
state 0(1 t) + I i)) that does not break the Z 2 symmetry. 
When 6 7 ^ 0 and a = c, the MPS is a symmetry breaking 
state of the form 141) + 17I'll) where U is Z 2 symmetry 
transformation. But when b ^ 0 and a c, what is the 
nature of the MPS? 

To answer such a question, we would like to study the 
symmetry structure of the environment matrix E. We 
find that, depends on which phase we are in, the symme¬ 
try structure of E will be very different. 

As we mentioned before, the environment matrix E 
and the associated scaling factor A is calculated via the 
iteration (or the self consistent condition) in Fig. 3. In 
general, there can be many environment matrices E that 
satisfy the self consistent condition. Here we choose those 
with largest absolute value of the scaling factor A. If 
there are many environment matrices with the degener¬ 
ate largest absolute value of the scaling factor, we then 
choose the environment matrices with minimal “entropy” 


5' = ^-Si In Si, ( 6 ) 

i 


where Si is the singular values of E. This will give us a 
set of environment matrices {E}. 

Next, we want to point out that environment matrix 
has only internal indices, so for E, the symmetry trans¬ 
formation (3) translates to: 

E^Ul-E- Ug. (7) 

If the M’s are invariant under the symmetry transforma¬ 
tion (3), then the set of environment matrices {E} will 
be invariant under the above transformation (7). 

If the action of transformation (7) is trivial on the set 
of environment matrices {E}, then the MPS does not 
have symmetry breaking. If the action is non-trivial {i.e. 
generate a permutation of the set {if}), then the MPS, in 
general, has a symmetry breaking; but this is not guar¬ 
anteed. 
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FIG. 6. The entanglement density matrix can be calcu¬ 
lated from the environment matrix. (The correct entangle¬ 
ment density matrix should be calculated from the total the 
environment tensor ^ ^ .) 
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FIG. 7. The off-diagonal term in environment matrix is 
plotted here as a function of a —c and b. We can see the delta- 
function-like behavior when a = c. Note the dip at a = c and 
6 = 0: this point corresponds to the four-fold degeneracy in 
the double-tensor, which is in the symmetric phase. 

The reason for the complication is that there are zero- 
measure possibilities that some internal bond degrees of 
freedom completely decouple from the physical degrees 
of freedom. To fix this problem, we may consider the 
entanglement density matrix - defined in 

Fig. 6 . We say two environment matrices are equivalent 
if they generate the same Pmima -Let us use 
{if}/ ^ to denote the equivalent class of the environment 
matrices. Then if the action (7) is non-trivial on {if}/ 
then the MPS has a symmetry breaking. 

With the above general discussion, we now go to our 
numerical results for internal dimension D — 2. For the 
symmetric phase, using the iteration method in Figure 
3 and after energy minimization, we obtain a final E = 

(^0 gives an invariant E under eqn (7). As 

a result, the total environment tensor has a pure tensor 
product form 

E*°* =E®E, ( 8 ) 

For the symmetry-breaking phase, after minimizing 
entropy according to eqn. ( 6 ), depending on the ini¬ 
tial values of E, the iteration method would give either 
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FIG. 8 . A plot of magnetization and the “fixed-point” order 
paramter as a function of h/J. When we increase the internal 
dimension, we get a more accurate phase transition point. 
In both cases, we get a first-order phase transition. When 
internal dimension D = 2, order parameter is just the off- 
diagonal term of the environment matrix (represented by blue 
dots). When D = 4, order parameter is the norm of the 2x2 
off-diagonal block in the environment matrix (represented by 
red circles). The black line curve is (erf) = ±(2 — 2h)'^'^®. 

if Si = (P or ifS 2 = ( ^ ^ 1 , which transforms 

\P Pj \-P P J 

into each other under eqn (7). Both of these correspond 
to environment matrix of the fixed-point wavefunction, 
as explained later in this section. If we construct the 
total environment tensor 

^ ^91 ^92 ^ ^92 ^ ( 9 ) 

then the Z 2 symmetry is restored, but now the the total 
environment tensor does not have a pure tensor product 
form. 

We are now in the position to answer the question 
raised at the begining of this section. We now know that 
symmetry-breaking phase is signatured by a non-zero off- 
diagonal term in the environment matrix. As shown in 
Fig. 7, if we plot this off-diagonal term as a function 
of a — c and 6 in M (see eqn (5)), then we see that the 
system is only in symmetry breaking state when a = c 
and 6 yf 0. So when 6 yf 0 and a ^ c, the state is in the 
symmetric phase. 

Here we’ve also plotted the magnetization as a func¬ 
tion of h/J in Fig. 8 . Note that in the graph, we get a 
first-order phase transition for both D = 2 and D = A. 
This is because we required our matrices M’s in the 
MPS to have the Z 2 symmetry (recall eqn. (5)). This 
symmetry requirement favors the symmetric phase, be¬ 
cause symmetry-breaking phase requires a = c, so M is 
block-diagonalized, reducing its effective internal dimen¬ 
sion. Thus the phase transition point is shifted leftwards, 
leading to a first-order transition. As we increase the in- 
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FIG. 9. A self-consistent condition that environment ma¬ 
trix must satisfy, where A® is a scaling factor. The part 
enclosed by the shaded area is the symmetry twisted donble- 
tensor 


ternal dimension, we expect the phase transition point to 
approach h/ J = 1 from the left. 

In Fig. 8, we’ve also plotted the order parameter as a 
function of h/ J. Note that these “order parameters” do 
not vary as we change hjJ. This is because the environ¬ 
ment matrix is obtained from enough iterations that it re¬ 
ally corresponds to the fully renormalized wavefunction. 
The order parameter obtained from the environment ma¬ 
trix then corresponds actually to the order parameter at 
the fixed point, thus is always the same until hitting the 
phase transition point. 


IV. DETECTING PHASES OF MPS WITHOUT 
KNOWING THE TRANSFORMATION 
PROPERTY OF THE MATRICES 


the matrices defined via the n**' power of double-tensor 
(where T and T® are viewed as matrices) 

Pah,a/3 = rh"^.aa, Pfh.a/3 = (P")h/3,aa- (10) 

If P and P® have different singular values in large n limit, 
then the corresponding MPS break the symmetry explic¬ 
itly. 

If I A® I = |A|, then the two environment matrices P® 
and P are related by the symmetry transformation (see 

( 3 )) 

E = ESUg, or E-^E0 = Ul. (11) 

In fact, we have 

P = UlEUg (12) 

If Ug forms a projective representation of the symmetry 
group G, then the corresponding MPS does not break the 
symmetry and has a non-trivial SPT order protected by 
the on-site symmetry. If Ug forms a ID representation 
of the symmetry group G, then the corresponding MPS 
does not break the symmetry and has a non-trivial SPT 
order protected by translation symmetry (and the on-site 
symmetry). 

Let us apply the above approach to a MPS state of 
spin-1 chian, where the matrix M*, I = x,y,z are given 
by the Pauli matrices: = aK The double-tensor is 

given by (see Fig. 3) 


In the above, we have assumed that the matrices in 
the MPS has the symmetry and studied how to use the 
symmetry breaking of the environment matrix to de¬ 
tect the spontaneous symmetry breaking in MPS. How¬ 
ever, in many calculations, such as the density-matrix- 
renormalization-group (DMRG) calculation, the result¬ 
ing matrices in the MPS do not have the symmetry in 
the symmetry breaking phase, and in general we do not 
even know how the matrices transform under the sym¬ 
metry transformation (since we do not know how the 
internal indices should transform under the symmetry). 
In this section, we will discuss how to detect the sponta¬ 
neous symmetry breaking in MPS, without knowing how 
the matrices in the MPS transforms under the symmetry 
transformation. 

Assume we have already obtained the matrices in the 
MPS. We first calculate the environment matrix E and 
the scaling factor A using Fig. 3. In general, the environ¬ 
ment matrix E is unique even in the symmetry break¬ 
ing state, since the matrices in the MPS obtained from 
DMRG in general already break the symmetry. Next, 
we insert the symmetry transformation g (see (3)) in the 
double-tensor to obtain a twisted double-tensor. The cor¬ 
responding twisted environment matrix is denoted as E^ 
and the twisted scaling factor as A® (see Fig. 9). 

If I A® I < |A|, then the corresponding MPS have a spon¬ 
taneous symmetry breaking. In fact, there is a more di¬ 
rect way to detect symmetry breaking. Let P and P® be 


E’bjS^aa 




(13) 


The action of the double-tensor Tbg^aa on the environ¬ 
ment matrix Ebp —>■ Tbp^aaEaa can be written in a matrix 
form 


E^^u^EuK (14) 

I 


We see that E = (the 2-by-2 identity matrix) is 

the non-degenerate environment matrix with A = 3. 

Now, let us show that the MPS has a x Z| symmetry 
where Z’f is generated by - the 180° spin rotation in 
S^^-direction and Z| is generated by Rz - the 180° spin 
rotation in S^-direction. Under the symmetry twists Rx 
and Rz, the corresponding double-tensors are 


rpRx _ -X X 

^bl3,a(x ^ba^j3ot 

rpRz _ _L 

^b0,aa '^ba^^a 


V V , z z 

- + ^baC^/3a 

y y z z 

^ba^^a. ^ba^^Oi' 


(15) 


The corresponding twisted environment matrices are 
given by 


= 2-i/V’", E^’‘ = 


(16) 


with = A^^ = —3. We see that 

Ur^ = Ur^ = (17) 
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FIG. 10. Applying imaginary time evolution in a layered 
structure. 


Since \X^- \ = |A^-| = |A| and 

E = EUr^ , E = UI^ EUr^ , (18) 

we found that the Z 2 x symmetry is not broken. We 
also see that Ur^, Ur^ generate a projective representa¬ 
tion of Z 2 X Z|. So the MPS is a SPT state protected 
by Zf X Z|. 


V. A TENSOR NETWORK APPROACH FOR 
ID MODEL 

In this section, we are going to use an infinite time- 
evolving block decimation (iTEBD) approach^'' to study 
ID models, such as the transverse Ising model (1). We 
are going study symmetry breaking by testing if |A®| = 
|A| or |A9| < |A|. 

A. The iTEBD method 

The iTEBD method is a tensor network version of 
the DMRG approach. The fundamental idea behind the 
iTEBD method is to use imaginary time evolution to get 
the ground state of a two-body Hamiltonian, and to use 
Singular Value Decomposition (SVD) to control the in¬ 
ternal dimensions. 

Consider any ID Hamiltonian with only nearest- 
neigbour interactions, we can always separate it into two 
parts, labeled by Ha and Hr- 

H = E 

i iGodd iGeven 

= Ha + Hb. (19) 

This way, either Ha or Hr would have no overlapping 
terms within itself. When the time step St is very tiny, 
we have: 


(20) 

We could then apply imaginary-time evolution layer by 
layer, as shown in Fig. 10. Now within each layer, time- 
evolution only operates on non-overlapping neighboring 
sites. Thus the entire problem reduces to a two-site prob¬ 
lem. 

The two-site time-evolution is done through Singular 
Value Decomposition, see Fig. 11. We first apply the 


■ j-s K 



Um’a.p 



m’ 

k’ 
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FIG. 11. Time-evolution on two sites. Step 1: apply time- 
evolution operator. Step 2: apply SVD and truncate the sin¬ 
gular matrix to only contain D largest singular values (S' de¬ 
notes the singular matrix after truncation). Step 3: seperate 
the singular values into the two sites. 


time-evolution operator (labeled by Wi^i+i) on two sites, 
resulting in a rank-4 tensor, Tm'a,k'y' 

Tm'a,k''i = Wm'k>,mkA^pBp~^- ( 21 ) 

Then we do SVD to split the rank-4 tensor: 

a,k'^ ~ Ujji' ( 22 ) 


Note that after applying SVD, the internal dimension 
has grown on the inner link. We could get back our 
original internal dimension by keeping only the D largest 
singular values of S'. We’ll call the truncated matrix 
S. Finally, we absorb diagnomal matrix S into the two 
on-site matrices, thus completing one step of evolution: 


4m _ TT 

^a/3 — IXir 


IdxjSp, Bp =Vki/d,j\ISi3. (23) 


This completes one step in time evolution. 

One improvement can be made on the above time- 
evolution step. Note that in the truncation process 
above, we implicitly assumed that all bond indices are 
equally important; however, we know that’s not the case. 
The “environment indices” a, 7 do not contribute equally, 
and their weights could be naturally included by us¬ 
ing the singular values Sa , S~f from the previous time- 
evolution step. (This is because in the previous step, a 
and 7 were inner link indices, and each index naturally 
carries a weight according to the previous step of SVD.) 

Thus the improved time-evolution step works as fol¬ 
lows: 1. First, we scale the “environment indices” using 
singular values Sa obtained from last step: 

A'Sp^^aA^p, (24) 

2. Then apply the time-evolution step described before. 

3. Lastly, we scale the “enviornment indices” back, by 
doing the following: 


^ Bf^ ^ (25) 


B. The iTEBD results for transverse Ising model 

After applying the imaginary-time evolution steps de¬ 
scribed above, we will eventually obtain the tensor that 
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FIG. 15. for the spin-1 model as a function of Jx 

with Jy = Jz = 1, Jzz = 0.4, calculated by the iTEBD method 
with D — 8. The “-I-” points are for symmetry twist g = Rx, 
“X” for g = Ry, and for g ^ Rz- When = 0, 

the corresponding symmetry g is not broken. We see that 
the phase near Jx = 1 has the full symmetry. The phase for 
smaller Jx has only the Ry symmetry. The phase for larger 
Jx has only the Rx symmetry. 


FIG. 12. The three up-down symmetric curves describe the 
magnetization of transverse Ising model (erf) as a function of 
h with J = 1. The “-I-” points are for internal dimension 
D = 1, for D — 2, and filled-box for D = 10. The other 
three curves describe The “x” points are for D — 1, 

for D = 2, and open-circle near (1,0) are for D — 10. 
The line curve is (erf) = ±(2 — 2/i)^^®. 



FIG. 13. The difference of the scaling factors of 

transverse Ising model as a function of h with J = 1. The 
“-I-” points are for D = 2, “x” for D = 4, “4^” for D = 8, 
and for D = 16. 



FIG. 14. The magnetization of transverse Ising model (af) 
as a function of h with J = 1. The “-I-” points are for D — 2, 
“x” for D — 4, for D = 8, for D = 16, and for 
D = 32. The line curve is (af) = ±(2 — 2/i)^/*. 



FIG. 16. The phase diagram for the spin-1 model in Jx-Jzz 
plane, with Jy = Jz = 1, calculated by the iTEBD method 
with D = 8. The points mark which = 0. The “-I-” 

points are for symmetry twist g = Rx, “x” for g = Ry, and 
for g — Rz. The green shaded area and the white area at 
the top have the full Rx, Ry, Rz symmetry. The gold shaded 
area has only the Ry symmetry, and the blue shaded area has 
only the Rx symmetry. 


describes the ground state wave function very well. The 
next issue is to identify the symmetry breaking order 
and/or SPT order in the ground state, using the method 
discussed before. 

For the transverse Ising model, Fig. 12 describes the 
calculated and the magnetization {a-}, using the 

iTEBD approach with various D. Fig. 13 and Fig. 14 
are the results near the transition point. The transition 
point is found to be he « 1.08 for D = 2, he ^ 1.0188 for 
D = 4, he^ 1.0101 for D = 8,he^ 1.0047 for D = 16, 
and he ~ 1.0015 for D = 32. The exact transition point is 
at he = 1. We see that “order parameter” works 

very well, in identify symmetry breaking transitions. 








































FIG. 17. A tensor-network state. All the physical sites are 
represented by dots, physical degrees of freedom by vertical 
lines, and internal degrees of freedom by in-plane links. 

C. The iTEBD calculation of ID model with 
symmetry-breaking and/or SPT orders 

In this section, we are going to use the iTEBD ap- 
praoch to study spin -1 model with x Z| symmetry: 

H = + jusm 

(26) 

The Z 2 is generated by 180° spin-rotation around 
the axis. The .Z| is generated by 180° spin-rotation 
Rz around the axis. The RxRz = Ry is the 180° 
spin-rotation around the Sy axis. 

We choose D = 8 and calculated the tensor for the 
ground state. To determine if describes a symmetry 
breaking state or not, it is not correct to directly test if 
has the symmetry or not. This is because even when 
is not invariant under any symmetry transformation 
of the form (3), can still describe a symmetric state. 

So to determine symmetry of the ground state, we 
instead calculated the quantity for symmetry 

twists g = Rx,Ry,Rz (see Fig. 15). We determined 
the phase diagram by examine where and which 
vanishes. The gold shaded area in Fig. 16 only has the 
Rx symmetry, since only = 0. The blue shaded 

area in Fig. 16 only has the Ry symmetry since only 
= 0 - Ttis green shaded area and the white area 

have = 0 for g = Rx,Ry,Rz and have the full 

Rx, Ry,Rz symmetry. In fact the green area is a phase 
with a non-trivial SPT order (the Haldane phase). 

VI. APPLICATION TO 2D MODEL WITH 
SYMMETRY BREAKING TRANSITION 

Now we want to generalize the above simple picture to 
2D. Consider the transverse Ising model on an infinite 
honeycomb lattice, with spins living on vertices. The 
Hamiltonian remains the same as equation (1). Since the 
ground state of a gapped system in 2D could be faithfully 
described by a tensor-network state,for a translation 
invariant system, we have (see Figure 17): 

I'J') = E E |{mj) (27) 

{mi} {a,/3,7} 


FIG. 18. The average energy with total environment tensor 
in 2D. Here, consists of four environment matrices, 
surrounding the two physical sites. 



FIG. 19. The self-consistent iteration process for environ¬ 
ment in 2D. 


where tensors M’s are again labeled by the physical de¬ 
grees of freedoms rrii, and tTr (tensor trace) contracts 
over all internal degrees of freedom on connected links 
labeled by a,f3 and 7 . Again, we want to do variational 
calculations with a simple picture involving the total en¬ 
vironment tensor which now consists of four envi¬ 

ronment matrices, see Figure 18. 

The key question now is how do we obtain a good envi¬ 
ronment matrix, as we did in ID (recall Figure 3)? Here 
we introduce a simple yet powerful iteration process: as¬ 
sume we have a three-fold rotational symmetry for ten¬ 
sor M, then the iteration needs two input matrices, and 
gives out only one output, see Figure 19. As before, af¬ 
ter enough numbers of iterations, we would reach a final 
stable “environment matrix”. 

It might be surprising, at first sight, why such a naive 
iteration process would give a reliable environment ma¬ 
trix. The key however is to realize that this iteration ac¬ 
tually gives an environment matrix for the infinite Bethe 
lattice, which is a very good first approximation for our 
honeycomb lattice (see Figure 20). As shown in the 
graph, the iteration process is actually equivalent to a 
self-consistent update for a large cluster of lattice points, 
and thus its legitimacy. 



FIG. 20. The iteration process on Bethe lattice. Note that 
this is a top-down view, and we have only shown one layer of 
tensor-network state. Starting from an environment surroud- 
ing C tensor, we could iterate to get environment for B, and 
then to A. The circled region is the region of interest. 
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FIG. 21. Energy as a function of h/J for the 2D transverse 
Ising model. As in the ID case, h- and J- terms are individu¬ 
ally plotted in the graph as well. We can see from the graph 
that our simulation shows a weak first-order phase transition, 
with phase transition point at h/J = 2.09. 



Just like in the ID case, here we would also like to 
require our tensor-product state to have a Z 2 symmetry 
corresponding to the spin up-down symmetry. Similar to 
the matrix-product state, on-site symmetry of the ground 
state also requires tensor M in the tensor-network state 
to transform in a special way: 




M m' _ ie, 

a./3,7 - e - 


Oc' ,0' , 7 ' 

(28) 


This is just a tensor generalization of condition (3). As 
before, m is the physical spin label, gmm' represents the 
on-site symmetry and acts in the spin space, Ug forms a 
projective representation of the symmetry group g, and 
is a unitary matrix acting on internal degrees of freedom 
labeled by a, /3 and 7 . 

For internal dimension of I? = 2, we can choose Ug = 

The most general symmetric tensor ^ 

satisfying eqn. (28) has 4 variational parameters and 
looks like the following: 

»i:r = {; )) =(‘ j). 




FIG. 22. Magnetization and order parameter as a function 
of h/J for the 2D transverse Ising model. When D — 2, 
the order parameter plotted (represented by red crosses) is 
p/ {s + t), which goes to zero in the symmetric phase. 


phase transition point occured at h/J = 2.09, which was 
within 2% error from Quantum Monte Carlo prediction 
of h/J = 2.13.'*® The typical runtime on a laptop was 
just a few seconds. If we increase the internal dimension 
to Z? = 4 and use a completely symmetric tensor with 
10 variational paramters, then we get a phase transition 
point at h/J = 2.12, within 1% error from the afore¬ 
mentioned Quantum Monte Carlo calculation. Here all 
variational parameters are real. 

Following what we did in iJ = 1, here we would also 
like to comment on the symmetry structure of the envi¬ 
ronment matrix E. Recall that E is obtained through 
iterations (or self consistent condition) in Fig. 19, which 
picks out the E with the largest absolute value of scal¬ 
ing factor A. One important difference/simplification in 
2D is that unlike in ID, in general, we do not have any 
degeneracies for E through the iteration equation (Fig. 
19), since the equation is non-linear. Thus in general E 
obtained is unique, and we do not need eqn. ( 6 ) to fix 
the basis. 


With the above discussion, we can go into the sym¬ 
metry structure for our environment matrix E (See Fig. 
22). For internal dimension being 2, using the iteration 
method mentioned above and after energy minimization, 

we have in the symmetric phase i? = M , which gives 


an invariant E under eqn. (7). As a result, the total en¬ 
vironment tensor is just the direct-product of them: 


Here we again assume M’s to be symmetric, because of 
rotational symmetry. 

With the above M tensors, numerical simulation could 
again be run on the 2D Ising model. Following what 
we did in ID, we vary h/J in equation (1) and mini¬ 
mize the energy for each value of h/J. By plotting the 
two energy terms, a phase diagram could also be ob¬ 
tained (see Fig. 21). For internal dimension D = 2, the 


E*°* = E®E®E®E. 


(30) 


As for the symmetry breaking phase, depending on the 


initial values of E, we get either E^^ = 


s p 
P i 


or E32 = 


s —p 
-P i 


which transforms into each other under eqn. 
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FIG. 23. Tensor-product state with spins living on links. 
Here physical sites are represented by dots. 

(7). If we construct the total environment tensor 

='^E3<g)ES<g,ES^E3, (31) 

g 

then the Z 2 symmetry is restored, but now again the total 
environment tensor does not have a pure tensor product 
form, as was the case in \D. 

VII. APPLICATION TO 2D MODEL WITH 
TOPOLOGICAL ORDER 

We now move on to the non-trivial example of Toric- 
Code model in a B-field, with spins living on links of an 
infinite 2D honeycomb lattice. The Hamiltonian is as 
follows: 

H = ( 32 ) 

V i^v p jGp k 

where cr’s are the usual Pauli matrices. We will first 
consider the phase diagram by fixing A —>■ 00 and varying 
h/B between [0, 00 ]. When h/B — 0, we have the original 
Toric-Code model, whose ground state is an equal-weight 
superposition of all closed loops of down-spins (in the 
background of up-spins). When h/B —>■ 00 , we have the 
spin-polarized state where all spins are pointing up. 

Later in this section, we will also consider the case 
when h/B —)■ —00 so the ground state is the fully packed 
loop state, which is an equal weight superposition of all 
loop configurations that are fully packed (every vertex 
has a loop passing through). The question of whether 
the fully packed loop state has topological order or not 
will then be explored. 

As in the previous example, we now try to use a tensor- 
network state to represent the ground state of the above 
Hamiltonian. Here since all spins live on links of the 
lattice, we will need two tensors T and M to represent 
our variational ground state (see Figure 23): 

= E E tTr[®„r„.; 3 ,^ M^/] |{m,}) (33) 

{mi} {a,/3,7,(5,T/} 

where v labels different vertices, I labels different links, 
a, (3, 7 , (5, rj label internal degrees of freedom, mi label 
physical degrees of freedom of link /, and tTr contracts 
over all connected internal indices. Note that due to the 
B term in the Hamiltonian 32, we will need to include an 



FIG. 24. Variational energy for Toric-Gode model in a B- 
field. 

entire plaquette in our variational calculation, as shown 
in Figure 24. 

Now we start by introducing our tensor ansatz in the 
simple case of internal dimension 2. In order to enforce 
the condition that A —>■ 00 and the rotational symmetry 
of the system, we need the following tensors T and M: 

f 1, 1/ a = /3 = 7 = 0; 

Ta,p,-y = \ X, else if a + l3 + ^ = 0 (mod 2); 

I 0 , otherwise; 

< = (J S) • < = (S ^ (3'*) 

where spin-up and spin-down’s are labeled by arrows. 
Note that when x = 1, it represents the regular Toric- 
Code ground state'^^, whereas when a; = 0, it represents 
the all-spins-up state. 

Before going into our variational calculations, we first 
note that our model in equation (32) could be mapped 
into a transverse Ising model by introducing a new pla¬ 
quette spin operator fip, where spins live on the plaque- 
ttes and p is the plaquette label.By doing the following 
mapping: HjGp ^ hp, MpMpo and consider only 

the A —)► 00 sector, our Hamiltonian reduces to: 

H = -Bj2h;-h ( 35 ) 

p <p,p'> 

which is the familiar transverse Ising model. Note that 
this Ising model is now on a 2D triangular lattice. 

With the above tensor network ansatz, we could run 
our variational scheme on the Toric-Code model. Just 
like in the Ising model cases, we vary h/B in eqn (32) 
(recall that we hold A —^ 00 ) and minimize the energy 
for each value of h/B. The environment tensor was cal¬ 
culated in the same way as before (See Fig. 19). The only 
difference here is that since the Hamiltonian (32) have a 
six-body interaction term, we have to include more sites 
into our mean-field calculation(See Figure 24). By plot¬ 
ting the two energy terms as a function of h/B, we get 
a phase diagram, which is plotted in Fig. 25. For in¬ 
ternal dimension D = 2, we got a phase transition point 
at B/h — 3.33, with an error of 30% to the Quantum 
Monte Carlo result of B/h = 4.768.'^^ This is not sur¬ 
prising as we only have one variational parameter. With 
internal dimension of 3 and only two variational param¬ 
eters, our result quickly improved to a phase transition 
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FIG. 25. Energy as a function of h/B for the 2D Toric Code 
model in a magnetic field. Here we only plotted the region 
when h/B > 0. We can see from the graph that our simu¬ 
lation shows a weak first-order phase transition, with phase 
transition point at h/ J = 0.3. 



FIG. 26. Internal Z 2 symmetry of the tensors T and M for 
our Toric-Code model. Here T has a (J 2 (g) cr^ (g) Uz symme¬ 
try, and M has a Uz (g) Uz symmetry, both of which square 
to identity. Note that this symmetry transformation could 
independently act on top and bottom layers of the tensor- 
network, thus the environment matrix transforms under a 
Z 2 X Z 2 group. 


point at B/h = 4.407, with an error of less than 8% to the 
Quantum Monte Carlo result. Note that the Quantum 
Monte Carlo value was obtained on the mapped equiva¬ 
lent model (see equation (35)) on a 2D triangular lattice. 

Now in order to understand the above result better 
and to further explore the case when h/B < 0, we need 
to understand the symmetry structure of both our tensor- 
product state and the environment tensor obtained. Note 
here that although the ground state doesn’t have a phys¬ 
ical Z 2 symmetry, the tensor ansatz T and M (34) still 
need to have an internal Z 2 symmetry, namely the “nec¬ 
essary symmetry condition” (See Figure 26): 

= (36) 

5 ' ,77' 

Here, the internal symmetry is represented by (7z®<Jz®(^z 
for tensor T, and CTz (g) Uz for tenor M, where (Tz is the 
Pauli matrix. Since both symmetry actions square to 
identity, we refer to the above internal symmetry as a Z 2 


symmetry. 

The physical reason for tensor ansatz to have the above 
“necessary symmetry condition” is that we want to make 
sure local variations of the tensors correspond to local 
perturbations of the Hamiltonian. Tensors that violates 
the above condition correspond to non-local perturbation 
in their Hamiltonian and thus can not be used to describe 
physical phase transitions"^®. 

It’s easy to check that tensors in equation (34) have the 
above symmetry. We would then like to ask, with the T 
and M tensors satisfying eqn (36) , what is the symmetry 
structure of the environment matrix? Note that unlike 
the Ising model, here the internal symmetry of the two 
layers of our tensor-network can act independently, as 
shown in Figure 26. Thus the environment matrix no 
longer transforms under eqn (7), but transforms under a 
Z 2 X Z 2 group: 


E-^Ul'E, or E-Ug, or Ul'E-Ug. (37) 


As in the Ising model, we expect that in different phases, 
the environment matrices E^s are either invariant under 
the above transformation, or undergoes a permutation. 

Our numerical result indeed shows the above feature. 
When internal dimension D = 2, we use the tensor ansatz 
in eqn (34) and iteration process (see Fig. 19) to get the 
environment matrix E. In the confined phase (including 

spin-polarized state), we obtain E = qV which is 


invariant under eqn. (37). As a result, the total environ¬ 
ment tensor is just the direct product of them: 


E*°^ = E®E®E®E®E®E. 


In the deconfined phase (including string-net state). 


ever, we have either E^^ = 


s 0 
0 t 


or ES^ = 


which transforms into each other under eqn. (37). 
could again construct a total environment tensor 


how- 



We 


Etot = (g) E3 (g) E3 (g) ES (g) ES (g) ES 

9 


that respects the Z 2 x Z 2 symmetry, but it does not have 
a pure tensor product form, as was the case for Ising 
model. 

In doing the above, we have really constructed a nu¬ 
merical way to detect topological orders. In the partic¬ 
ular case above, Z 2 topological order is signatured by a 
“symmetry breaking” in the environment matrix, which 
breaks the original Z 2 x Z 2 symmetry of E (see eqn (37)) 
down to Z 2 (see eqn (7)). 

With this realization, a natural question to ask is: if 
we now consider negative magnetic field with h/B —>■ 
— 00 , will the fully packed loop state has Z 2 topological 
order? To answer this question, let us first write down 
the ground state wave function of the fully packed loop 
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FIG. 27. The above graph shows a “string crystal” state, 
where blue lines represent spin-downs forming vertical strings, 
and red lines represent spin-ups forming the background. 

state in tensor form: 

f 0, if a = P = j = 0; 

Ta,i3,-y = t if a + P + "/ = 0 (mod 2 ); 

I 0 , otherwise; 



FIG. 28. The self-consistent iteration process for environ¬ 
ment matrix when the tensors do not have rotational symme¬ 
try. We introduce three different environment matrices, E^, 
and E'^. One cycle of iteration consists of three steps: 
1. Input E^ and E^ to update E^; 2. Input E^ and E^ 
to update E^-, 3. Input E^ and E^ to update E^. After 
iterating for enough number of steps, all three self-consistent 
equation will be simultaneously satished. 


= 

0 , 7 ] 



Mj = 

0 , 7 ] 



(38) 


Note the difference between this and eqn (34): here we 
require loops to cover each vertex, so Tg^o.o = 0 - 

Now to see whether this state has Z 2 topological order 
or not, all we need to do is to calculate its environment 
matrix through iteration (see Fig. 19). Depending on the 

initial condition, we obtain either = 

or Q ^ 0^62^ ’ again transforms into 

each other under eqn (37). This means that we are still 
in the deconfined phase, and packed loop state has Z 2 
topological order. 

One may worry that the simple test above would fail to 
differentiate the “string crystal” state (see Fig. 27) where 
string configuration is stationary, from the fully packed 
loop state where the string configurations are fluctuating. 
This worry turns out to be unnecessary through careful 
study below. 

Consider a “string crystal” state, with vertical strings 
formed by down-spins (shown in Fig. 27). We would like 
to study the symmetry structure of the environment ma¬ 
trix for this state. Note that unlike the previous tensor- 
network ansatz in eqn (34) and (38), here the tensors no 
longer have three-fold rotational symmetry. 

Our method could be easily generalized to non- 
rotationally symmetric tensors by introducing a three- 
step iteration process, shown in Fig. 28. (Recall that 
this is different from the symmetric iteration process in 
Fig. 19.) In this three-step iteration process, we intro¬ 
duce three different environment matrices, which then 
iterate in a cyclic fashion. Now, any physical quantities 
could again be calculated by sandwiching the operator in 
between two layers of tensor network states, surrounded 
by three different types of environment matrices shown in 


/0.38 0 \ 

0 0.62j 



FIG. 29. Expectation value of an operator with non- 
rotational-symmetric tensor networks. Note here that there 
are three different types of environment matrices. 


Fig. 29, and all of our previous analysis follows. (Again 
compare this with Fig. 18, where there was only one type 
of environment matrix.) 

With the above three-step iteration process, the envi¬ 
ronment matrix of our “string crystal” state could then 
be obtained as follows: 

=(J S)'-®'’==(D!) ■ 

Note that the above E^ and E^ do not really break the 
Z 2 X Z 2 symmetry shown in eqn (37). This is because 
the iteration process for E^ and E'^ are both linear, so 
an overall minus sign does not affect the iteration result. 
Thus E^ and E^ could only be determined up to a sign, 
which is a fictitious gauge degree of freedom and have no 
physical meaning. Thus E^ —E^ does not correspond 
to breaking the Z 2 x Z 2 symmetry, and “string crystal” 
state indeed does not possess Z 2 topological order. 


VIII. CONCLUSIONS 

In this paper, we proposed a new signature for phase 
transitions between tensor-network states using the en- 
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vironment matrices. Different phases are distinctively 
labeled by different symmetry structures in the environ¬ 
ment matrices. Thus through carefully studying differ¬ 
ent symmetry structures of the environment matrix, we 
could identify different phases and obtain the detailed 
phase boundaries for both symmetry-breaking transitions 
and topological phase transitions. This greatly helps us 
in identifying topological orders or SPT orders from a 
generic tensor-network state. 

The environment matrix is obtained through a very 
simple iteration process using the tensor-network state, 
in both ID and 2D. This iteration process provides a 
self-consistent environment matrix that summarizes the 
contributions from far away sites, and is like a “mean- 
field” theory for tensor-networks. In the same line of 
thinking, the environment matrix serves like an “order 
parameter”. What’s special about this “mean-field” the¬ 
ory is that it’s suitable for studying long-range entangled 
states, and is thus suitable for tackling topological phase 
transitions. 

In ID, we demonstrated that this new signature could 
be easily combined with existing numerical methods like 
DMRG or iTEBD to identify SPT phases. We first ob¬ 
tain the ground state in a matrix-product form by ap¬ 
plying these ID numerical methods. Then we calculate 
the environment matrix, either through direct iteration 
process or through a twisted iteration process (where the 


symmetry transformation g is sandwiched in between the 
double tensor in the iteration process). By simply com¬ 
paring the scaling factors in the two iteration process, we 
could identify which SPT phase we are in, thus provid¬ 
ing an easy way to identify SPT orders directly from a 
matrix-product state. 

In 2D, the iteration process gives a very efficient way 
of calculating variational energies, which in turn leads 
to a simple numerical methods in obtaining gound state 
wave function by minimizing the energy. If we require the 
ground state tensors to have the proper on-site symmetry, 
iteration process could give us environment matrices that 
have drastically different symmetry structures, labeling 
different (topological) phases. Note that the on-site sym¬ 
metry doesn’t have to be a physical symmetry— internal 
gauge symmetry is also valid. 

The above numerical method is very general and could 
be easily applied to many interesting systems in higher 
dimension including 3D systems. This will open new 
doors in numerical study of higher dimensional systems. 
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